SYDE 556/750: Simulating Neurobiological Systems

Terry Stewart

Spatial Representation

  • A common feature that researchers find in the activity of real neurons is that there's some sort of spatial structure in how the neurons are organized
    • Topographic maps
  • Example: retinotopic maps (the layout of neurons maps onto some spatial organization from the retina):

  • Example: auditory cortex organized by frequency sensitivity

  • Example: superior colliculus for eye movement

  • How do we fit this into the NEF?

Spatial structure and encoders

  • One aspect is nice and easy
  • Instead of randomly choosing encoders, organize them in a nice way
    • the location of the neuron relates to its encoder
  • Nothing else changes
  • There's even an option in the interactive mode interface to automatically do this (even while a model is running)
  • Does this solution handle everything?
    • What's missing?

Consider vision

  • Consider the visual example
    • The simplest approach would be to say each neuron has some preferred point in space it responds to
    • 2-dimensional representation $x = <h, v>$
    • each neuron gets an encoder $e$ based on its position
    • $a = G[\alpha e \cdot x + J^{bias}]$
  • That's not quite going to work
    • a point $x$ farther away from the origin, but in the same direction as $e$ will make this neuron fire faster
    • need to do a change of variables and add an extra dimension
    • $x = <h, v, \sqrt{1 - h^2 -v^2}>$
  • Now I can give a point as input and the neurons in the correct area of the visual cortex will fire
  • What's missing here?

Complex representations

  • But when you see things, you get more than one input at a time
  • What happens if I see a dot at $<h_0, v_0>$ and another dot at $<h_1, v_1>$ at the same time?
  • Is that the same as seeing one dot at $<h_0 + h_1, v_0 + v_1>$?
    • No
  • We see them as two distinct points
  • Same thing with auditory (multiple frequencies at once, plus other sound features)
  • Maybe even the same thing with the superior colliculus and eye movements (although that's less clear)

  • So what should this representation be?

    • Not just 2-dimensional
    • Need to represent a whole image

Grid-based representation

  • We know how to represent vectors, so why not just break the visual area up into pixels, and have one dimension per pixel?
    • So for a 3x3 image of a square, the $x$ value might be $<1,1,1,1,0,1,1,1,1>$
    • Can do brightness (and maybe even colour) this way too
  • Now choose encoders based on location
  • Simplest approach: each encoder looks like $<0,0,0,1,0,0,0,0,0>$
    • This is equivalent to nine separate ensembles
    • What are the limitations of this approach?
    • How can we improve this?

Changing how we choose encoders

  • Perfectly aligned encoders: $<0,0,0,1,0,0,0,0,0>$
  • Randomly chosen encoders: $<0.3, -0.1, 0.02, 0.1, 0.7, -0.2, -0.3, 0.001, -0.2>$
  • Can also go inbetween: $<0.25, 0, 0, 0.5, 0.25, 0, 0.25, 0, 0>$
  • What is that?
    • A neuron that responds mostly to dots in one location, but also a little bit to dots nearby
  • Why would we do this?
    • What computational advantage do we get?

Computation with sparse local encoders

  • With perfectly aligned encoders, we can only do linear functions of $x$
  • With randomly chosen encoders, we can do any function of $x$, but the accuracy decreases very quickly as we get up to higher-order terms (since there are many pairs of variables)
  • With encoders that are local (i.e. all the non-zero terms are clustered together) and sparse (i.e. fairly few non-zero terms), we can compute fairly complex function, as long as they are local (i.e. they're dependent on things that are spatially nearby)
    • So we'll be worse at things that require us to relate two points far away in space, but better at things where the points are nearby
    • This seems like a reasonable tradeoff

Does the brain do this?

  • Looks to be the case
  • Very common way to look at neurons in the visual field
  • Find the stimulus a neuron most responds to

  • Usually modelled as a gabor filter (a 2-D sine wave times a Gaussian):

  • Pretty much any vision machine learning algorithm gives things like this

  • So this could be a good way to generate encoders for visual stimuli
    • And maybe even other stimuli too

Implementation

  • Let's start with making gabors to use as our encoders

In [6]:
import numpy

def gabor(width, height, lambd, theta, psi, sigma, gamma, x_offset, y_offset):
    x = numpy.linspace(-1, 1, width)
    y = numpy.linspace(-1, 1, height)
    X, Y = numpy.meshgrid(x, y)
    X = X - x_offset
    Y = Y + y_offset

    cosTheta = numpy.cos(theta)
    sinTheta = numpy.sin(theta)
    xTheta = X * cosTheta  + Y * sinTheta
    yTheta = -X * sinTheta + Y * cosTheta
    e = numpy.exp( -(xTheta**2 + yTheta**2 * gamma**2) / (2 * sigma**2) )
    cos = numpy.cos(2 * numpy.pi * xTheta / lambd + psi)
    return e * cos

g = gabor(500, 500, lambd=0.6, 
                    theta=numpy.pi/4, 
                    psi=numpy.pi/2, 
                    sigma=0.2, 
                    gamma=0.7,
                    x_offset=0,
                    y_offset=0)
import pylab
pylab.imshow(g, extent=(-1, 1, -1, 1), interpolation='none', vmin=-1, vmax=1, cmap='gray')
pylab.show()


  • Now make a bunch of them

In [7]:
width = 50
height = 50

N = 100

def make_random_gabor(width, height):
    return gabor(width, height, 
                  lambd=random.uniform(0.3, 0.8),
                  theta=random.uniform(0, 2*numpy.pi),
                  psi=random.uniform(0, 2*numpy.pi),
                  sigma=random.uniform(0.2, 0.5),
                  gamma=random.uniform(0.4, 0.8),
                  x_offset=random.uniform(-1, 1),
                  y_offset=random.uniform(-1, 1))

encoders = [make_random_gabor(width, height) for i in range(N)]
import pylab
pylab.figure(figsize=(10,8))
for i in range(N):
    w = i%12
    h = i/12
    pylab.imshow(encoders[i], extent=(w, w+0.95, h, h+0.95), interpolation='none',
                 vmin=-1, vmax=1, cmap='gray')
    pylab.xticks([])
    pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, N/12))
    
pylab.show()


  • now let's make some neurons and find some decoders
  • to find decoders, we need to choose $x$ values to sample over
    • what should we choose?

In [8]:
S = 100
width = 50
height = 50

x = numpy.random.normal(size=(width*height, S))

import pylab
pylab.figure(figsize=(10,8))
for i in range(S):
    w = i%12
    h = i/12
    img = x[:,i]
    img.shape = width, height
    pylab.imshow(img, extent=(w, w+0.95, h, h+0.95), interpolation='none',
                 vmin=-1, vmax=1, cmap='gray')
    pylab.xticks([])
    pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, S/12))
    
pylab.show()


  • Why might that be a bad idea?
  • Other options?

In [9]:
S = 100

width = 50
height = 50

def normalize(v):
    return v/numpy.linalg.norm(v)


C = 4
x = [normalize(numpy.sum([make_random_gabor(width, height) for i in range(C)],axis=0).flatten()) for j in range(S)]
x = numpy.array(x).T

import pylab
pylab.figure(figsize=(10,8))
for i in range(S):
    w = i%12
    h = i/12
    img = x[:,i]
    img.shape = height, width
    pylab.imshow(img, extent=(w, w+0.95, h, h+0.95), interpolation='none',
                 vmin=-0.1, vmax=0.1, cmap='gray')
    pylab.xticks([])
    pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, S/12))
    
pylab.show()



In [10]:
import syde556

S = 100
x = []
for i in range(S):
    sx, A = syde556.generate_signal(T=width, dt=1.0, rms=1, limit=0.05, seed=i)
    sy, A = syde556.generate_signal(T=height, dt=1.0, rms=1, limit=0.05, seed=i+2000)
    SX, SY = numpy.meshgrid(sx, sy)
    img = SX*SY
    x.append(normalize(img.flatten()))
x = numpy.array(x).T


import pylab
pylab.figure(figsize=(10,8))
for i in range(S):
    w = i%12
    h = i/12
    img = x[:,i]
    img.shape = height, width
    pylab.imshow(img, extent=(w, w+0.95, h, h+0.95), interpolation='none',
                 vmin=-0.1, vmax=0.1, cmap='gray')
    pylab.xticks([])
    pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, S/12))
    
pylab.show()


  • Now let's build the neurons and see how good they are

In [11]:
import syde556
reload(syde556)

N = 100
width = 50
height = 50

encoders = [make_random_gabor(width, height).flatten() for i in range(N)]

vision = syde556.Ensemble(neurons=N, dimensions=width*height, encoders=encoders)

decoder = vision.compute_decoder(x, x, noise=0.1)

A, xhat = vision.simulate_rate(x, decoder)

In [12]:
import pylab
pylab.figure(figsize=(10,12))
for i in range(S):
    w = i%6
    h = i/6
    img = x[:,i]
    img.shape = height, width
    pylab.imshow(img, extent=(2*w, 2*w+0.9, h, h+0.9), interpolation='none',
                 vmin=-0.1, vmax=0.1, cmap='gray')
    img = xhat[:,i]
    img.shape = height, width
    pylab.imshow(img, extent=(2*w+0.95, 2*w+1.85, h, h+0.9), interpolation='none',
                 vmin=-0.1, vmax=0.1, cmap='gray')
    pylab.xticks([])
    pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, S/6))
    
pylab.show()


Pixel size

  • What size of pixel should we use?
  • How does that affect things?

In [17]:
width = 50
height = 50

# make samples
S = 100
x = []
for i in range(S):
    sx, A = syde556.generate_signal(T=width, dt=1.0, rms=1, limit=2.5/width, seed=i)
    sy, A = syde556.generate_signal(T=height, dt=1.0, rms=1, limit=2.5/height, seed=i+2000)
    SX, SY = numpy.meshgrid(sx, sy)
    img = SX*SY
    x.append(normalize(img.flatten()))
x = numpy.array(x).T
print x.shape

#make neurons
N = 500

encoders = [make_random_gabor(width, height).flatten() for i in range(N)]

vision = syde556.Ensemble(neurons=N, dimensions=width*height, encoders=encoders)

decoder = vision.compute_decoder(x, x, noise=0.1)

A, xhat = vision.simulate_rate(x, decoder)


scale = 5/np.sqrt(width*height)
#plot results
import pylab
pylab.figure(figsize=(10,12))
for i in range(S):
    w = i%6
    h = i/6
    img = x[:,i]
    img.shape = height, width
    pylab.imshow(img, extent=(2*w, 2*w+0.9, h, h+0.9), interpolation='none',
                 vmin=-scale, vmax=scale, cmap='gray')
    img = xhat[:,i]
    img.shape = height, width
    pylab.imshow(img, extent=(2*w+0.95, 2*w+1.85, h, h+0.9), interpolation='none',
                 vmin=-scale, vmax=scale, cmap='gray')
    pylab.xticks([])
    pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, S/6))
    
pylab.show()


(2500L, 100L)
  • Why is this the case?
  • What would happen if we went up to an infinite number of pixels?
  • How could we implement this?

Representing Continuous Spaces

  • It seems unnatural to pick some particular number of pixels
    • especially since actual vision doesn't really have anything like pixels
  • So what can we do instead?
  • Let $x$ be a continuous function, rather than a vector
  • What does that mean?
    • Instead of representing $x = <0, 0.5, 1, 0.5, 0, -0.5, -1, -0.5>$
    • we might have $x = sin(v)$ where $v \in [0, 2\pi)$
  • What else would have to change?
    • $e$ would also be a function (over the same range of $v$ values)
    • instead of $x \cdot e$, we would do $\int x(v) e(v) dv$
    • normalize so that $\int e(v) e(v) dv=1$
  • Problems:
    • How do we do that integral?
    • Would have to do it analytically
    • Seems restrictive
    • If we do it analytically and approximate the integral with some finite $dv$ then we're doing exactly the same thing as if we had have had some large (but finite) number of pixels
  • How else can we represent a function?

Representing a function

  • Here's another way to represent a function:
    • $f(v) = \sum_i c_i b_i(v)$
    • where $b_i(v)$ is the set of basis functions
  • Example: the polynomials:
    • $b_i(v) = v^i$
    • so I might have $f(v)= v^3 + 3 v^2 + 3v +1$
    • so the coefficients $c$ are $<1, 3, 3, 1>$
    • This means we can represent $f(v)= v^3 + 3 v^2 + 3v +1$ using the vector $<1, 3, 3, 1>$
    • And we know how to represent vectors with the NEF
  • Is that all we need?
  • Still have to be able to compute $\int x(v) e(v) dv$

  • $\int x(v) e(v) dv$

  • $\int (\sum_i c_{x,i} b_i(v)) (\sum_i c_{e,i} b_i(v)) dv$
  • $\int (\sum_i c_{x,i} c_{e,i} b_i(v)^2)dv + (\mbox{ ... all the cross terms ... })$
  • so those cross terms are ugly
    • can we get rid of them?
  • Let's assume $b_i(v)$ are orthogonal

    • $\int b_i(v) b_j(v) dt$ is zero for $i \neq j$
    • all the cross terms disappear
  • $\int (\sum_i c_{x,i} c_{e,i} b_i(v)^2)dv$

  • now rearrange, pulling the sum out
  • $\sum c_{x,i} c_{e,i} \int b_i(v)^2 dv$
  • what if we make the basis functions have unit length?
    • $\int b_i(v)^2 dv = 1$
  • so $\int x(v) e(v) dv$ becomes $\sum c_{x,i} c_{e,i}$, which is just $c_x \cdot c_e$
    • so our dot product operation on the functions just becomes exactly the dot product of the vector of coefficients
    • so we won't have to change any of our code at all!
  • the only works assuming the basis vectors $b_i(v)$ are:
    • orthogonal: $\int b_i(v) b_j(v) dt=0$ for $i \neq j$
    • normalized: $\int b_i(v)^2 dv = 1$
  • this is called an orthonormal basis

An othonormal basis

  • Let's go back to the polynomial example
  • $b_i(v) = v^i$
  • Are these orthogonal?
    • $\int b_i(v) b_j(v) dt=0$ for $i \neq j$
    • $\int v^i v^j dt=0$ for $i \neq j$
    • $\int v^{i+j} dt=0$ for $i \neq j$
    • let's pick a range of $v \in [-1,1]$
    • ${1 \over {i+j+1}} [(1)^{i+j+1} - (-1)^{i+j+1}]$
    • that'll only be zero if $i+j+1$ is even....

In [134]:
import numpy
v = numpy.linspace(-1, 1, 1000)
b_i = v**1
b_j = v**2

import pylab
pylab.plot(v, b_i, label='$b_i$')
pylab.plot(v, b_j, label='$b_j$')
pylab.plot(v, b_i*b_j, label='$b_i b_j$')
pylab.legend(loc='best')
xlabel('$v$')
pylab.show()
dv = 2.0/len(v)

print 'sum:', numpy.sum(b_i*b_j*dv)


sum: 3.90312782095e-18
  • So what could we use instead?
  • There is an orthogonal basis set that's fairly similar to the polynomials
    • We've even seen them before
    • The Legendgre polynomials

In [135]:
for i in range(5):
    plot(v, numpy.polynomial.legendre.legval(v, numpy.eye(5)[i]))
xlabel('$v$')

show()


  • let's make sure they're orthogonal

In [136]:
import numpy
v = numpy.linspace(-1, 1, 1000)
i = 1
j = 2
b_i = numpy.polynomial.legendre.legval(v, numpy.eye(i+1)[i])
b_j = numpy.polynomial.legendre.legval(v, numpy.eye(j+1)[j])

import pylab
pylab.plot(v, b_i, label='$b_i$', linewidth=3)
pylab.plot(v, b_j, label='$b_j$', linewidth=3)
pylab.plot(v, b_i*b_j, label='$b_i b_j$')
pylab.legend(loc='best')
xlabel('$v$')
pylab.show()
dv = 2.0/len(v)
print 'sum:', numpy.sum(b_i*b_j*dv)


sum: 9.28077059648e-17
  • Are they normalized?

In [137]:
v = numpy.linspace(-1, 1, 1000)
dv = 2.0/len(v)

length = []
for i in range(10):
    b_i = numpy.polynomial.legendre.legval(v, numpy.eye(i+1)[i])
    length.append(sum(b_i*b_i)*dv)

plot(length)  
xlabel('$i$')
show()


  • Nope.
  • The area between -1 and 1 for $b_i$ is $2/(2i+1)$
  • So we rescale

In [138]:
v = numpy.linspace(-1, 1, 1000)
dv = 2.0/len(v)

length = []
for i in range(10):
    b_i = numpy.polynomial.legendre.legval(v, numpy.eye(i+1)[i])
    b_i = b_i * numpy.sqrt((2*i + 1)/2.0)
    length.append(sum(b_i*b_i)*dv)

plot(length)  
xlabel('$i$')
ylim(0, 2)
show()


  • What do these scaled basis functions look like?

In [294]:
for i in range(5):
    plot(v, numpy.polynomial.legendre.legval(v, numpy.eye(i+1)[i])* numpy.sqrt((2*i + 1)/2.0))
xlabel('$v$')    
show()


  • So, let's say I want a neural population that can represent polynomials
  • Here are a set of randomly generated $x$ values showing example functions to represent

In [295]:
D = 5
S = 20
x_poly = numpy.random.randn(D, S)/numpy.sqrt(D)

v = numpy.linspace(-1, 1, 1000)

for i in range(S):
    plot(v, numpy.polynomial.polynomial.polyval(v, x_poly[:,i]))
xlabel('$v$')    
show()


  • Now we convert those into our basis space

In [296]:
x = []
for i in range(S):
    x.append(numpy.polynomial.legendre.poly2leg(x_poly[:,i])/numpy.sqrt((2*i + 1)/2.0))
x = numpy.array(x).T

v = numpy.linspace(-1, 1, 1000)

for i in range(S):
    plot(v, numpy.polynomial.legendre.legval(v, x[:,i])*numpy.sqrt((2*i + 1)/2.0))
xlabel('$v$')
show()


  • Now let's make a neural population to represent these
  • It's just a perfectly normal 5-dimensional population
    • We still pick encoders randomly

In [297]:
N = 50

ensemble = syde556.Ensemble(neurons=N, dimensions=5)

decoder = ensemble.compute_decoder(x, x, noise=0.1)

A, xhat = ensemble.simulate_rate(x, decoder)

figure(figsize=(10,4))
subplot(1, 2, 1)
for i in range(S):
    plot(v, numpy.polynomial.legendre.legval(v, x[:,i])*numpy.sqrt((2*i + 1)/2.0))
ylim(-3, 3)
xlabel('$v$')
title('original')
subplot(1, 2, 2)
for i in range(S):
    plot(v, numpy.polynomial.legendre.legval(v, xhat[:,i])*numpy.sqrt((2*i + 1)/2.0))
ylim(-3, 3)
xlabel('$v$')
title('decoded')
show()


  • We can do representation
  • Can we do computation?
  • Same way as before
    • All we have to do is define what function we want to approximate
  • Let's start with computing the area under the curve

In [298]:
def area(x):
    fv = numpy.polynomial.legendre.legval(v, x)*numpy.sqrt((2*i + 1)/2.0)
    a = sum(fv)*dv
    return [a]

fx = numpy.array([area(xx) for xx in x.T]).T
    
decoder = ensemble.compute_decoder(x, fx, noise=0.1)

A, fxhat = ensemble.simulate_rate(x, decoder)

for i in range(10):
    print x[:,i], fx[:,i], fxhat[:,i]


[-0.45249468  0.06689621  0.65862321 -0.22781363  0.04283414] [-3.99012265] [-3.692147]
[-0.13508849 -0.45085403 -0.29704591 -0.14351425 -0.04635415] [-1.19610606] [-0.99389806]
[ 0.03477385  0.59657275 -0.18919625  0.00530379 -0.06779654] [ 0.30484099] [ 0.29829337]
[ 0.15631005 -0.07939436 -0.10217138 -0.09528619 -0.02525089] [ 1.37936599] [ 1.31180309]
[-0.45944895  0.02760009 -0.27678184  0.08308695 -0.0547291 ] [-4.06067516] [-3.88563593]
[ 0.25147024  0.0522454   0.02486886  0.11510253  0.03265593] [ 2.22143421] [ 1.61950481]
[ 0.03977521 -0.06903214  0.03796262  0.09403689  0.02458351] [ 0.35183856] [ 0.24418011]
[-0.34478934  0.26520947 -0.15590749  0.04113701 -0.00300367] [-3.04650196] [-2.60464404]
[-0.08437981  0.28442028  0.19314713  0.05949096 -0.01813703] [-0.74367552] [-0.55534047]
[-0.01543183 -0.11761898  0.06478489  0.07507399 -0.00886165] [-0.13579599] [-0.09870697]
  • What about finding the maximum?

In [300]:
def maximum(x):
    fv = numpy.polynomial.legendre.legval(v, x)*numpy.sqrt((2*i + 1)/2.0)
    index = numpy.argmax(fv)
    return [v[index]]

fx = numpy.array([maximum(xx) for xx in x.T]).T
    
decoder = ensemble.compute_decoder(x, fx, noise=0.01)

A, fxhat = ensemble.simulate_rate(x, decoder)

for i in range(20):
    print x[:,i], fx[:,i], fxhat[:,i]


[-0.45249468  0.06689621  0.65862321 -0.22781363  0.04283414] [-1.] [-1.03426409]
[-0.13508849 -0.45085403 -0.29704591 -0.14351425 -0.04635415] [-0.94394394] [-0.94790649]
[ 0.03477385  0.59657275 -0.18919625  0.00530379 -0.06779654] [ 0.78178178] [ 0.80341805]
[ 0.15631005 -0.07939436 -0.10217138 -0.09528619 -0.02525089] [ 0.21721722] [ 0.03715061]
[-0.45944895  0.02760009 -0.27678184  0.08308695 -0.0547291 ] [-0.17317317] [-0.14223553]
[ 0.25147024  0.0522454   0.02486886  0.11510253  0.03265593] [ 1.] [ 0.98869827]
[ 0.03977521 -0.06903214  0.03796262  0.09403689  0.02458351] [ 1.] [ 0.47082208]
[-0.34478934  0.26520947 -0.15590749  0.04113701 -0.00300367] [ 1.] [ 1.01930901]
[-0.08437981  0.28442028  0.19314713  0.05949096 -0.01813703] [ 1.] [ 0.98959533]
[-0.01543183 -0.11761898  0.06478489  0.07507399 -0.00886165] [-0.7997998] [-0.36437875]
[-0.147665    0.049168   -0.06538418 -0.00584206 -0.02632643] [ 0.47347347] [ 0.32563004]
[ 0.02441503 -0.1699608   0.05092713 -0.02245458 -0.00107661] [-1.] [-0.90374256]
[ 0.14030728  0.06652304  0.04019052  0.00865063  0.0332155 ] [ 1.] [ 0.9777674]
[-0.04854667  0.19268431  0.12686373 -0.00477267 -0.00366835] [ 1.] [ 0.95878839]
[ 0.14449928  0.04155529  0.23186591  0.01521891  0.03536285] [ 1.] [ 0.94144684]
[-0.04972338 -0.09333049  0.195611   -0.03597476  0.02794989] [-1.] [-0.93327958]
[ 0.0649419  -0.02901023  0.03064615 -0.03862747 -0.00013708] [-1.] [-0.87645674]
[ 0.11407232 -0.05903258  0.02590225 -0.0141109   0.01488021] [-1.] [-0.69223594]
[ 0.07999296  0.02683351  0.0958427   0.03251879  0.02420397] [ 1.] [ 0.71611585]
[ 0.13212116  0.19791227  0.16534859  0.02255834  0.03551369] [ 1.] [ 1.22114539]
  • We can also compute functions as output
    • $f(x) = -x$
    • But, since $x$ is a function, this is really:
    • $f(x(v))=-x(v)$
  • This one's really easy, since we can do it right on the basis functions

In [301]:
def negative(x):
    return -x

fx = numpy.array([negative(xx) for xx in x.T]).T
    
decoder = ensemble.compute_decoder(x, fx, noise=0.01)

A, fxhat = ensemble.simulate_rate(x, decoder)

figure(figsize=(10,4))
subplot(1, 2, 1)
for i in range(S):
    plot(v, numpy.polynomial.legendre.legval(v, x[:,i])*numpy.sqrt((2*i + 1)/2.0))
xlabel('$v$')
ylim(-3, 3)
title('original')
subplot(1, 2, 2)
for i in range(S):
    plot(v, numpy.polynomial.legendre.legval(v, fxhat[:,i])*numpy.sqrt((2*i + 1)/2.0))
xlabel('$v$')
ylim(-3, 3)
title('decoded')
show()


  • So it's all the same as normal NEF
  • What functions $x$ we're good at representing and what functions $f(x)$ we're good at computing depends on the basis functions
  • How do we find these for things other than the polynomials?

Making your own othonormal basis functions

  • Let's go back to the vision example
  • We need a set of basis functions that will be good at encoding $x$ (images) and $e$ (the gabor-shaped tuning curves for the neurons)
  • So, given a bunch of functions, we need to find an orthonormal basis space
  • How do we do that?
  • make a big set of the functions, evaluated at some $dv$
  • do singular value decomposition
  • let's start with a simple set of tuning curves: gaussian bumps

In [187]:
N = 200

v = numpy.linspace(-1, 1, 1000)
def make_gaussian(centre, sigma):
    return numpy.exp(-(v-centre)**2/(2*sigma**2))

encoders = numpy.array([make_gaussian(centre=random.uniform(-1, 1), 
                          sigma=0.1) for i in range(N)])
plot(v, encoders[:20].T)
xlabel('$v$')
show()


  • Now do SVD

In [192]:
U, S, V = numpy.linalg.svd(encoders.T)

plot(v, U[:,:10])
xlabel('$v$')
show()


  • Let's make sure they're orthogonal

In [198]:
import numpy
v = numpy.linspace(-1, 1, 1000)
i = 0
j = 1
b_i = U[:,i]
b_j = U[:,j]

import pylab
pylab.plot(v, b_i, label='$b_i$', linewidth=3)
pylab.plot(v, b_j, label='$b_j$', linewidth=3)
pylab.plot(v, b_i*b_j, label='$b_i b_j$')
pylab.legend(loc='best')
xlabel('$v$')
pylab.show()
dv = 2.0/len(v)
print 'sum:', numpy.sum(b_i*b_j*dv)


sum: -2.50635984098e-20
  • Are they normalized?

In [209]:
v = numpy.linspace(-1, 1, 1000)
dv = 2.0/len(v)

length = []
for i in range(10):
    b_i = U[:,i]
    length.append(sum(b_i*b_i)*dv)

plot(length)  
xlabel('$i$')
ylim(0, 0.004)
show()


  • Not quite
  • Missing the $dv$ factor

In [224]:
v = numpy.linspace(-1, 1, 1000)
dv = 2.0/len(v)

basis = U/(math.sqrt(dv))

length = []
for i in range(10):
    b_i = basis[:,i]
    length.append(sum(b_i*b_i)*dv)

plot(length)  
xlabel('$i$')
ylim(0, 2)
show()


  • Here's what they look like after being scaled

In [228]:
plot(v, basis[:, :10])
xlabel('$v$')
show()


  • How many dimensions should we use?
  • The SVD tells us that too
  • The singular values indicate how important each basis function is for accurately representing the desired functions

In [223]:
plot(S)
show()


  • So we'll pick the first 20 or so to be our basis
  • In order to use these, we need to convert our encoders $e$ and our samples $x$ to be coefficients in this basis space

In [288]:
bases = 20
figure(figsize=(10,4))
subplot(1, 2, 1)
plot(v, encoders.T[:,:20])
ylim(0,1)
title('original')

subplot(1, 2, 2)
basis = U[:,:bases]/(math.sqrt(dv))
e = numpy.dot(encoders, basis)*dv
plot(v, numpy.dot(basis, e.T)[:,:20])
ylim(0,1)
title('bases=%d'%bases)
show()


  • That gives us 20-dimensional encoders that represent these gaussian functions well
  • We can do the same thing to some random $x$ functions we want these neurons to represent

In [289]:
S = 100
import syde556
x = numpy.array([syde556.generate_signal(2.0, 2.0/1000, 
                  rms=0.5, limit=2, seed=i)[0] for i in range(S)])

plot(v, x[:5].T)
xlabel('$v$')
show()


  • Here's how well that basis space represents those functions

In [290]:
figure(figsize=(10,4))
subplot(1, 2, 1)
plot(v, x.T[:,:20])
title('original')

subplot(1, 2, 2)
basis = U[:,:bases]/(math.sqrt(dv))
x_vector = numpy.dot(x, basis)*dv
plot(v, numpy.dot(basis, x_vector.T)[:,:20])
title('bases=%d'%bases)
show()


  • Now we can do standard NEF
    • The coefficients of the encoder functions $e$ are our encoders
    • The coefficients of the sample functions $x$ are our samples

In [332]:
N = 200
bases = 20
S = 200 

encoders = numpy.array([make_gaussian(centre=random.uniform(-1, 1), 
                          sigma=0.1) for i in range(N)])

U, Sing, V = numpy.linalg.svd(encoders.T)

basis = U[:,:bases]/(math.sqrt(dv))
e = numpy.dot(encoders, basis)*dv

x_func = numpy.array([syde556.generate_signal(2.0, 2.0/1000, 
                  rms=0.5, limit=2, seed=i)[0] for i in range(S)])

x = numpy.dot(x_func, basis).T*dv

ensemble = syde556.Ensemble(neurons=N, dimensions=bases, encoders=e)

decoder = ensemble.compute_decoder(x, x, noise=0.1)

A, xhat = ensemble.simulate_rate(x, decoder)


figure(figsize=(10,4))
subplot(1, 2, 1)
plot(v, x_func.T[:,:20])
ylim(-4, 4)

title('original')

subplot(1, 2, 2)
basis = U[:,:bases]/(math.sqrt(dv))
xhat_func = numpy.dot(basis, xhat)
plot(v, xhat_func[:,:20])
title('bases=%d'%bases)
ylim(-4, 4)
show()


  • We can also compute functions on this representation

In [347]:
N = 200
bases = 20
S = 200 

encoders = numpy.array([make_gaussian(centre=random.uniform(-1, 1), 
                          sigma=0.1) for i in range(N)])

U, Sing, V = numpy.linalg.svd(encoders.T)

basis = U[:,:bases]/(math.sqrt(dv))
e = numpy.dot(encoders, basis)*dv

x_func = numpy.array([syde556.generate_signal(2.0, 2.0/1000, 
                  rms=0.5, limit=2, seed=i)[0] for i in range(S)])

x = numpy.dot(x_func, basis).T*dv

ensemble = syde556.Ensemble(neurons=N, dimensions=bases, encoders=e)

def my_function(x):
    value = numpy.dot(x, basis.T)
    value = -value
    #value = numpy.hstack([value[100:], value[:100]])
    return numpy.dot(value, basis)*dv
    
    
fx = numpy.array([my_function(xx) for xx in x.T]).T    


decoder = ensemble.compute_decoder(x, fx, noise=0.1)

A, fxhat = ensemble.simulate_rate(x, decoder)


figure(figsize=(10,4))
subplot(1, 2, 1)
plot(v, x_func.T[:,:20])
ylim(-4, 4)

title('original')

subplot(1, 2, 2)
basis = U[:,:bases]/(math.sqrt(dv))
fxhat_func = numpy.dot(basis, fxhat)
plot(v, fxhat_func[:,:20])
title('bases=%d'%bases)
ylim(-4, 4)
show()



In [ ]: